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Abstract 

We  propose  a  new  metamodeling  method  to  characterize 
the  output  (response)  random  process  of  a  dynamic  system 
with  random  parameters,  excited  by  input  random  processes. 
The  metamodel  can  be  then  used  to  efficiently  estimate  the 
time-dependent  reliability  of  a  dynamic  system  using  analytical 
or  simulation-based  methods.  The  metamodel  is  constructed 
by  decomposing  the  input  random  processes  using  principal 
components  or  wavelets  and  then  using  a  few  simulations  to 
estimate  the  distributions  of  the  decomposition  coefficients.  A 
similar  decomposition  is  also  performed  on  the  output  random 
process.  A  kriging  model  is  then  established  between  the  input 
and  output  decomposition  coefficients  and  subsequently  used 
to  quantify  the  output  random  process  corresponding  to  a 
realization  of  the  input  random  parameters  and  random 
processes.  What  distinguishes  our  approach  from  others  in 
metamodeling  is  that  the  system  input  is  not  deterministic  but 
random.  The  quantified  output  random  process  is  finally  used 
to  estimate  the  time-dependent  reliability  or  probability  of 
failure  of  the  dynamic  system  using  the  total  probability 
theorem.  The  proposed  method  is  illustrated  with  a  numerical 
example. 

1.  Introduction 

The  response  of  dynamic  systems  under  uncertainty  is 
described  by  a  random  process.  The  input  commonly  consists 
of  a  combination  of  random  variables  and  random  processes. 
A  time-dependent  reliability  analysis  is  therefore,  needed  to 
calculate  the  probability  that  the  system  will  perform  its 
intended  function  successfully  for  a  specified  time.  It  is 
therefore,  related  to  product  functionality  over  time. 

Reliability  is  an  important  engineering  requirement  for 
consistently  delivering  acceptable  product  performance 
through  time.  As  time  progresses,  the  product  may  fail  due  to 
time-dependent  operating  conditions  and  material  properties, 
component  degradation,  etc.  In  this  research,  we  use  time- 
dependent  reliability  concepts  associated  with  the  first-passage 
of  non-repairable  systems. 

The  time-dependent  probability  of  failure,  or  cumulative 
probability  of  failure  [1 , 2],  is  defined  as 


Pf  (o,  t)  =  p{ar  e  [o,  t]  :  g(p,  z(f),  t)  >  St  }  (1 ) 

where  g  is  a  random  process  that  maps  the  input  random 
variables  (3  and  random  processes  z(r)  to  a  response  of 
interest  and  St  is  a  given  threshold  value.  The  random 

process  G(t)=g(^z(l)j)-St  can  be  viewed  as  a  collection 

of  random  variables  at  different  time  instances  f.  Since  we 
consider  a  first  excursion  failure  problem  in  Equation  (1),  the 

failure  domain  is  defined  as  F  =  j  ijiax^g(p,z(f),f)  >  St  j. 

The  time-dependent  probability  of  failure  of  Equation  (1) 
can  be  calculated  exactly  as 

Pf  (o,  t)  =  1  -  (1  -  Pf  (0))  exp  j-  u(r)df|  (2) 

where  Tj.- (0)  is  the  instantaneous  probability  of  failure  at  the 

Pit  <  Tc  <  t  +  dr  I  Tf  >  t) 

initial  time  and  Mt)  =  lim  -  is  the 

dr— >0  dr 

failure  rate  with  Tf  denoting  the  time  to  failure.  In  the 

commonly  used  out-crossing  rate  approach,  the  failure  rate  is 
approximated  by  the  up-crossing  rate 

p[g(p,Z(0,r)-S  <0ng(p,Z(r  +  Ar),r  +  Ar)-S  >  oj 

v+(r)=  lim  - - - - - 

Ar  ->  0,  Ar 

Ar  >  0 

(3) 

under  the  assumptions  that  the  probability  of  having  two  or 
more  out-crossings  in  [r,  r  +  Ar]  is  negligible,  and  the  out- 

crossings  in  [r,r  +  Ar]  are  statistically  independent  of  the 

previous  out-crossings  in  [o,  r]. 

Monte  Carlo  simulation  (MCS)  can  accurately  estimate  the 
probability  of  failure  in  Equation  (1)  but  it  is  computationally 
prohibitive  for  dynamic  systems  with  a  low  failure  probability.  To 
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address  the  computational  issue  of  MCS,  analytical  methods 
have  been  developed  based  on  the  out-crossing  rate  approach 
which  was  first  introduced  by  Rice  [3]  followed  by  extensive 
studies  [2,  4-6]  under  the  assumption  that  the  out-crossings 
are  statistically  independent  and  Poisson  distributed.  The  PHI2 
method  [2]  uses  two  successive  time-invariant  analyses  based 
on  FORM,  and  the  binomial  cumulative  distribution  to  calculate 
the  probability  of  the  joint  event  in  Equation  (3).  A  Monte-Carlo 
based  set  theory  approach  has  been  also  proposed  [7,  8]  using 
a  similar  approach  with  the  PHI2  method.  Because,  analytical 
studies  (e.g.,  [9-11])  have  shown  that  the  PHI2  based  approach 
lacks  sufficient  accuracy  for  vibratory  systems,  analytical 
approaches  have  been  proposed  to  accurately  estimate  the 
time-dependent  probability  of  failure  considering  non¬ 
monotonic  behavior  [9,  12,  13]. 

The  limited  accuracy  of  the  out-crossing  rate  approach 
has  been  recently  improved  by  considering  the  correlations 
between  the  limit-state  function  at  two  time  instants  [11].  The 

method  estimates  the  up-crossing  rate  v+(t)  by  solving  an 
integral  equation  involving  v+(t)  and  v++(t,  fj),  the  joint 
probability  of  up-crossings  in  times  t  and  rl  [14]. 

Among  the  simulation-based  methods,  a  MCS  approach 
was  proposed  in  [15]  to  estimate  the  time-dependent  failure 
rate  over  the  product  lifecycle  and  its  efficiency  was  improved 
using  an  importance  sampling  method  with  a  decorrelation 
length  [16]  in  order  to  reduce  the  high  dimensionality  of  the 
problem.  Subset  simulation  [17,  18]  has  been  recently 
developed  as  an  efficient  simulation  method  for  computing 
small  failure  probabilities  for  general  reliability  problems.  Its 
efficiency  comes  from  introducing  appropriate  intermediate 
failure  sub-domains  to  express  the  low  probability  of  failure  as 
a  product  of  larger  conditional  failure  probabilities  which  are 
estimated  with  much  less  computational  effort.  Because  it  is 
very  challenging  to  generate  samples  in  the  conditional 
spaces,  subset  simulation  with  Markov  Chain  Monte  Carlo 
(SS/MCMC)  [19]  and  subset  simulation  with  splitting  (SS/S) 
[20,  21]  methods  have  been  introduced. 

In  this  paper,  we  describe  a  new  metamodeling  method  to 
characterize  the  output  random  process  of  a  dynamic  system 
with  random  parameters,  excited  by  input  random  processes 
which  can  be  used  to  efficiently  estimate  the  time-dependent 
probability  of  failure  using  either  analytical  or  simulation-based 
methods.  The  metamodel  is  constructed  by  decomposing  the 
input  random  processes  using  principal  components  or 
wavelets  and  then  using  a  few  simulations  to  estimate  the 
distributions  of  the  decomposition  coefficients.  A  similar 
decomposition  is  also  performed  on  the  output  random 
process.  A  kriging  model  is  then  established  between  the  input 
and  output  decomposition  coefficients  and  subsequently  used 
to  quantify  the  output  random  process  corresponding  to  a 
realization  of  the  input  random  parameters  and  random 
processes.  The  proposed  approach  is  new  because  it 
considers  the  system  input  as  random  and  not  deterministic  as 
is  usually  the  case. 

The  paper  is  arranged  as  follows.  Section  2  describes  the 
proposed  metamodeling  methodology  and  Section  3  uses  a 
mathematical  example  to  demonstrate  its  applicability.  Finally, 
Section  4  provides  a  brief  summary  and  concludes. 


2.  PROPOSED  METAMODELING 
APPROACH 


Assume  that  the  simulation  model  has  two  types  of  inputs: 
random  time-dependent  inputs  represented  by  the  random 
processes  Z(r)and  random  time-independent  inputs 
represented  by  random  variables  (3.  Assume  that  the  output  is 
also  time-dependent  and  is  represented  by  the  random 
processes  Y(t).  We  decompose  each  process  zfrjusing  a 
well-defined  basis  of  functions  B(r)  such  as  wavelets  or 
principal  components.  For  principal  components,  consider  D 
realizations  of  the  random  process  Zfrjdenoted  by 

Zj( 1 1...,  zD(t).  An  eigenvalue  decomposition  of  their 
covariance  matrix  defines  the  principal  components.  If  we 
retain  a  small  number  rof  dominant  eigenvalues  and  (0  is  a 
vector  whose  rows  are  the  corresponding  eigenvectors,  each 
z  j(t\j  =  l,...,ocan  be  expressed  as 

Zj(t)=j.aiBu(t)  =  a  ;Bj  (t)  where a  ;-  =  .  A 

similar  principal  decomposition  can  be  considered  for  the 
output  process  Y(t)  so  that  each  output  realization  (sample 

function)  Y:(  t)  j  =  l, ....  D  can  be  expressed  as 

Yj(t)=isiB2i(t)  =  djB2  (t)  using  basis  functions  B,(oand 
T 

5  j  =  Yj(r)B2(t)  ,  if  we  retain  a  small  number  s  of  dominant 
eigenvalues  of  the  output  covariance  matrix. 

To  simplify  the  description  of  the  proposed  method,  we 
assume  that  we  have  only  one  input  random  variable  p,  one 
input  random  process  Z(t),  and  one  output  random  process 
Y(t)  which  is  represented  with  only  one  principal  component. 
The  input  vector  can  be  therefore,  expressed  in  terms  of  the 
vector  (a,  P)  where  a  =  {ay,  a2,  ar )  and  the  output  is 

expressed  in  terms  of^.  We  use  the  principal  component 

decomposition  in  this  paper  although  other  bases  may  be 
chosen  (e.g.  wavelets). 


What  distinguishes  our  approach  from  others  in 
metamodeling  is  that  the  input  vector  (a,  P)  is  not  deterministic 
but  random.  The  distributions  of  all  components  of  (a,  P)  are 
estimated  using  a  small  number  of  computer  runs.  If  we 
assume  that  all  a! s  are  independent  and  normally  distributed 
as  A/(yu, ,  cr,2)  where  /U,  is  the  mean  of  a,  over  all  performed 
computer  runs  and  a,  is  the  standard  error  of  cr,  over  all 
computer  runs,  and  also  assume  that  p  is  a  normal  random 
variable  with  mean  yUp  and  standard  error  op,  independent  of 
the  ci’s,  the  joint  probability  density  function  between  a  and  p  is 
given  by 


/«*,/?) 


1  I"  1  (P-Hp)2 

- <=ex p  - 2 - 

a p'iln  2  ap 


(4) 
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We  run  the  mathematical  model  for  D  input  vectors  (a(1), 
P(1)),  ...  ,  (a(D),  p(D))  and  obtain  the  corresponding  time- 

dependent  outputs  whose  first  principal  component  is  (5(1),..., 

S([\  Denote  the  vectors  a  =(a<1' . a(D|),  §=(P(1) . P(D))  and  5 

=  (5(1),.--,  5(D)).  Conditioned  on  (a,  |3),  the  vector  5  has  a 
multivariate  Gaussian  distribution  with  mean  r/ 1  and 
2 

covariance  matrix  z  C  where  C  is  a  Gaussian  correlation 
matrix  whose  elements  depend  on  the  distances  between 
sampled  inputs  (a,  P).  The  parameters  of  this  distribution  are 
estimated  by  maximum  likelihood. 

Let  (oiq,/?q)  be  a  new  input  vector  and  50  the 

corresponding  output  first  principal  component.  A  kriging 
metamodel  is  generated  between  the  input  (a,  |3)  and  the 
output  5  to  obtain  50.  Details  about  kriging  metamodeling  can 
be  found,  for  example,  in  [22-24],  The  kriging  Gaussian 

prediction  distribution  f(S()  I  a.  P)  for  <50,  is  given  by 


f(S0  I  a,  P)  =  .  exp 

P 

0  2/  T  -1  \ 

where  the  prediction  variance  is  va  p  =  r  \1-CqC  Cq  j 

and  the  prediction  mean  is  p  =  //  +  Cq  c  ^d-ni).  Here 

Cq  is  the  vector  of  correlations  between  the  new  input  and  the 
sampled  inputs.  The  joint  distribution  of  <?Q,a,p  is 


2v 


'_L  o 
0  ro  “  wa,p 
a,p 


(5) 


/ (<?0,  a,  p)  =  / (<?0  I  a,  P)/ (a,  P)  oc 

0  ^ 
a,P  ) 


D 

I!  exp 

7=1 

The  marginal  distribution  of  £q  at  the  new  input  vector 
(aQ,/^)  is  given  by  the  integral 

f(S0)  =  \f(S0,a,p)dadp.  (7) 

Note  that  the  typically  small  number  of  runs  performed  in 
an  engineering  problem  is  used  only  to  estimate  the 

distributions /(a,  p)  and  f(SQ  la,p)of  Equations  (4)  and  (5), 
respectively.  Once  these  distributions  are  estimated,  a  larger 
number  N  of  statistical  simulations  from  the  distribution  /(5q) 

of  Equation  (7)  can  be  obtained.  For  each  simulated  <Jq  we 
create  a  simulated  output  function  Y0(t)  =5qB^U)  ■ 


'  9  h  2  9  2 

2  ay  2  afj 


3U) 


(6) 


-  exp 


2m. 


a,p 


2v. 


a,p 


We  now  provide  some  details  on  the  evaluation  of  integral 
in  Equation  (7).  Consider  R  sampled  sets  of  D  input  vectors 

(«,  P)j,...,  («,  P)^ .  The  integral  is  approximated  using  the 

Monte  Carlo  method,  up  to  a  proportionality  constant,  as 

1  R  .  1  R  . 

-Tf(S0Aa,p)k)  =  -Tf(S0\(a,P)k)f((a,P)k).  (8) 

R  K— i  R  K—l 

While  /((a,  p)^ )  can  be  computed  easily  for  any  number  R  of 
sets,  f(SQ  I  (a,  P)^ )  will  require  running  the  mathematical 

model  R*D  times;  i.e.,  for  each  set  [(a,  p)^ ,  y]  of  D  input 

vectors.  This  is  computationally  very  demanding.  Instead,  we 
use  a  jackknife  method  (e.g.  [25]),  which  belongs  to  the 
general  class  of  “leave-one-out”  methods  for  cross-validation. 

(-k)  (-k) 

Specifically,  we  choose  R=D  and  («,  P),  ;=  («  ,p  ), 
meaning  that  we  replace  (a,  P),  by  the  set  of  D  original  input 

vectors  a,p  from  which  we  leave  the  input  vector  out,  for  k 

=  1 , . . . ,  D.  The  advantage  of  this  method  is  that  we  already  have 

(- k )  (- k ) 

data  from  math  model  runs  at  these  (a  ,p  )  vectors, 
and  we  do  not  need  any  additional  math  model  runs. 

The  total  probability  theorem  is  used  to  compute  the  time- 
dependent  probability  of  failure.  We  draw  M  new  random 

inputs  (aQ,/?Q)from  the  input  distribution  /(a,/?) of  Equation 
(4)  and  for  each  of  them  simulate  N  output  random  functions 
Fq(H by  sampling  /(Jq)  of  Equation  (7)  N  times.  For  each  of 

the  M  random  inputs  (ciq,  4)).  a  simple  Monte  Carlo 

simulation  is  used  to  calculate  a  probability  of  failure  P^1  (o,  t) 

using  the  N  output  functions  Yq(i).  It  can  be  easily  shown 
using  the  total  probability  theorem  that  the  probability  of  failure 
is  the  average  among  the  Py  (o,  t)  probabilities;  i.e., 

(  \  1  n  Mr  \ 

Pf(p,t)= — ZPf  \0 ,t).  Therefore,  the  time-dependent 

probability  of  failure  is  the  percentage  of  the  M*N  response 
functions  YQ(t)  that  exceed  a  specified  failure  threshold  T  in 

the  interval  [0,  f\.  This  probability  considers  the  randomness  of 
the  inputs  Zfrjand  p. 

3.  EXAMPLE 

Consider  the  following  mathematical  model 
y(/)  =  2exp(z(r))sin(2^»)  (9) 

where  Z(f)  is  a  random  process  generated  from  the  function  2- 
t,  to  which  correlated  noise  has  been  added.  The  correlated 
noise  is  represented  by  random  draws  from  a  Gaussian 
distribution  with  a  zero  mean  vector  and  Gaussian  correlation 
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for  the  temporal  dimension  of  smoothness  parameter  25.  The 
random  parameter  p  is  normally  distributed  as  A/(2,  0.05). 

A  sample  of  D  =  30  computer  runs  is  obtained  (Figure  1). 
Figure  la  shows  the  sample  of  30  random  input  functions  Z(f) 
for  t  e  [o,  2]  and  Figure  1b  shows  their  principal  component 
reconstruction  using  r  =  4  principal  components  which  account 
for  approximately  90%  of  the  variance.  Figure  1c  shows  the 
corresponding  30  random  output  functions  Y(f)  and  Figure  Id 
shows  their  principal  component  reconstruction  using  only  the 
first  principal  component  which  accounts  for  about  66%  of  the 
variability.  Visual  inspection  shows  that  the  time-dependent 
output  reconstruction  captures  the  major  features  of  actual 
time-dependent  output,  even  with  only  the  first  principal 
component. 

The  random  input  vector  is  expressed  with  5  components: 
the  4  coefficients  from  the  principal  component  decomposition 
of  the  input  function  Z(f)  denoted  by  a  and  the  p.  The  time- 
dependent  output  is  expressed  by  a  scalar  which  is  the 
dominant  principal  component  coefficient  5.  The  four 
components  of  a  are  assumed  independent  and  normal,  with 
estimated  means  equal  to  zero  and  standard  deviations  equal 
to  0.9062,  0.8387,  0.6854,  and  0.4947,  respectively.  These 
estimated  values  are  obtained  from  the  sample  of  30  computer 
runs.  The  estimated  mean  and  standard  deviation  of  p  are 
1.9997  and  0.0531  respectively,  which  are  in  agreement  with 
the  theoretical  values  of  2  and  0.05.  The  scatter  plots  of  Figure 
2  verify  the  assumption  that  the  four  components  of  a  are 
independent. 

We  have  developed  a  metamodel  using  the  proposed 
approach  using  M  =  100  new  random  inputs  drawn  from  the 
input  distribution  and  N  =  50  simulations  for  each  of  these  new 
random  inputs.  Therefore,  we  have  a  total  of  5,000  simulations 

TQfrjshown  in  green  in  the  left  panel  of  Figure  3.  We  observe 

that  the  variability  of  these  curves  appears  to  be  similar  with 
the  variability  of  the  30  actual  math  model  output  curves 
( Figure  1c)  or  their  principal  components  reconstructions 
( Figure  Id).  Flowever,  their  number  is  5,000  (much  greater 
than  30),  which  allows  for  a  more  accurate  evaluation  of  the 
probability  of  failure.  The  right  panel  of  Figure  3  shows  the 
time-dependent  probability  of  exceeding  the  thresholds  T  of  4, 
6,  7.5,  1 0  or  1 5.  We  observe  as  expected,  that  as  the  threshold 
is  increased  the  time-dependent  probability  of  failure  reaches 
its  maximum  value  of  1  at  a  later  time.  The  accuracy  of  the 
calculated  probabilities  of  failure  can  be  easily  verified  using 
Monte  Carlo  simulation. 
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Figure  1b.  Principal  component  reconstruction  using  r=  4 


Figure  1c.  Corresponding  30  random  output  functions  Y(t) 


Figure  la.  A  sample  of  30  random  input  functions  Z(f) 
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Figure  Id.  Principal  component  reconstruction 
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Figure  2.  Scatter  plots  of  the  four  components  of  vector  a 


Figure  3.  Metamodel  simulations  (left  panel)  and  time-dependent  probability  for  different  thresholds  (right  panel) 
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4.  Summary/Conclusions 

We  described  in  this  paper  a  new  metamodeling  approach 
to  characterize  the  output  random  process  of  a  dynamic  (or 
time-dependent  in  general)  system  with  random  input 
parameters,  excited  by  an  input  random  process.  The 
methodology  utilizes  a  decomposition  of  the  input  and  output 
random  processes  using  principal  components  or  wavelets  and 
kriging.  The  latter  establishes  a  non-parametric  interpolation 
between  the  input  and  output  decomposition  coefficients  which 
is  then  used  to  quantify  the  output  random  process.  The 
metamodel  is  used  to  efficiently  estimate  the  time-dependent 
probability  of  failure  using  either  analytical  or  simulation-based 
methods.  This  capability  is  very  practical  and  desirable  for 
large-scale,  linear  or  non-linear,  dynamic  systems  where  the 
calculation  of  response  is  time  intensive  and/or  expensive  to 
compute  or  measure.  The  proposed  approach  considers  the 
system  input  as  random  and  not  deterministic  as  is  usually  the 
case.  We  used  a  simple  mathematical  example  to  demonstrate 
the  steps  of  the  approach  and  demonstrate  its  practicality. 
Future  work  will  generalize  the  method  for  non-Gaussian 
random  variables  and  processes. 
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